Ghost introgression facilitates genomic divergence of a sympatric cryptic lineage in Cycas revoluta

Abstract A cryptic lineage is a genetically diverged but morphologically unrecognized variant of a known species. Clarifying cryptic lineage evolution is essential for quantifying species diversity. In sympatric cryptic lineage divergence compared with allopatric divergence, the forces of divergent selection and mating patterns override geographical isolation. Introgression, by supplying preadapted or neutral standing genetic variations, can promote sympatric cryptic lineage divergence via selection. However, most studies concentrated on extant species introgression, ignoring the genetic legacy of introgression from extinct or unsampled lineages (“ghost introgression”). Cycads are an ideal plant for studying the influence of ghost introgression because of their common interspecific gene flow and past high extinction rate. Here, we utilized reference‐based ddRADseq to clarify the role of ghost introgression in the evolution of a previously identified sympatric cryptic lineage in Cycas revoluta. After re‐evaluating the evolutionary independency of cryptic lineages, the group‐wise diverged single‐nucleotide polymorphisms among sympatric and allopatric lineages were compared and functionally annotated. Next, we employed an approximate Bayesian computation method for hypothesis testing to clarify the cryptic lineage evolution and ghost introgression effect. SNPs with the genomic signatures of ghost introgression were further annotated. Our results reconfirmed the evolutionary independency of cryptic lineage among C. revoluta and demonstrated that ghost introgression to the noncryptic lineage facilitated their divergence. Gene function related to heat stress and disease resistance implied ecological adaptation of the main extant populations of C. revoluta.


| INTRODUC TI ON
A cryptic lineage is a genetically distinct group of descendants whose superficial morphologies are indistinguishable from those of a known species. Recognition of this hidden diversity and its ecological role helps us to appreciate the real biodiversity of the world better and formulate workable conservation policies for reducing biodiversity loss arising from ecosystem stresses like climate change (Gomez-Corrales & Prada, 2020;Hodkinson et al., 2011). Because of the morphogenetic discordance of cryptic lineages, we are interested in how they evolved and how the lineage entities persist genetically (Michalski & Durka, 2015;Monro & Mayo, 2022;Struck et al., 2018).
Compared with allopatric cryptic lineages, which may be mainly influenced by geographical isolation, sympatric cryptic lineages with overlapping population distributions may be caused by divergent selection as explained by Titus et al. (2019), or assortative mating, or allopatric divergence following secondary contact, which was put forward by Binks et al. (2021). Natural selection dominates in the first two mechanisms and will drive adaptive divergence of specific genes, causing genomic islands (Foote, 2018). These genomic islands can have pleiotropic effects ("magic traits") connected with nonrandom mating that act as barrier loci to promote speciation. Other loci linked to genomic islands could hitchhike to produce genomewide divergence. These selection-driven processes will become efficient under high genetic variations, which can be facilitated by introgression.
Introgression from a third lineage may broaden ecological opportunities ("adaptive introgression") and confer standing genetic variations, advancing population divergence into new environments (Oziolor et al., 2019;Richards & Martin, 2017) or even facilitating divergence of the introgressed lineage from its sister lineage, making sympatric speciation easier (Richards et al., 2019). A mounting number of studies that focused on introgression between present species, however, ignored the genetic legacy left by extinct (archaic) or unsampled lineages, also known as ghost introgressions (Beerli, 2004). However, we maintain that consideration of ghost introgression is a requisite for correctly inferring the mechanisms of macro-and micro-evolutionary processes (Hahn & Nakhleh, 2016;Herrera-Alsina et al., 2022;Tricou et al., 2022aTricou et al., , 2022b. By transferring standing genetic variations to introgressed lineages, ghost introgression is another force for accelerating adaptation and speciation (Ottenburghs, 2020). Although ghost introgression is well known for facilitating human colonization and adaptation (Gokcumen, 2020;Racimo et al., 2015Racimo et al., , 2017, there has been much less focus on its application to wild plants. This is the first study to explore the genomic divergence of cycads from the perspective of ghost introgression, and Cycas revoluta is an outstanding system for investigating the influence of ghost introgression on the formation of sympatric cryptic diversity. Cycas revoluta is mainly distributed along the Ryukyu-Taiwan Archipelago, and a previous study first found a sympatric cryptic lineage called TaiB, restricted to Taiwan (Figure 1a) with unknown origin and functional divergence . The main coastal habitats of C. revoluta are ecologically different from those of its sister species, Cycas panzhihuaensis, from the inland Sichuan basin of China (Liu et al., 2018;Xiao & Moller, 2015), within the area where Cycas originated . The disjunct distribution of C. revoluta with long phylogenetic branches and morphological autapomorphy imply a high extinction-to-speciation rate during the evolution of C. revoluta (Condamine et al., 2015;Liu, Lindstrom, et al., 2022;Walters & Osborne, 2004). The occurrence of multiple extinction events together with frequent hybridization among Cycas (Norstog & Nicholls, 1997) makes it worth studying the influence of ghost introgression on sympatric cryptic lineage evolution in C. revoluta.
Leveraging the published C. panzhihuaensis genome (Liu, Wang, et al., 2022), we realigned the ddRADseq sequences from Chang et al. (2022) (Figure 1b, Table S1). They treated Cycas taitungensis, endemic to Taiwan, as a synonym for C. revoluta based on overlapping diagnostic traits, paraphyletic relationships, and larger intrathan interspecific genetic variations. However, a diverged cryptic lineage restricted to Taiwan was also discovered. According to their genetic and geographical distribution, we separated populations of C. revoluta into three groups: Ryu, TaiA, and the cryptic lineage, TaiB (Figure 1a). Group Ryu comprises the Ryukyu and Fujian populations (a small population on the China coast, genetically admixed with the Ryukyu populations); group TaiA is genetically mixed with Ryu, but geographically isolated in Taiwan; and group TaiB is a sympatric cryptic lineage with TaiA in Taiwan, which cannot be distinguished morphologically, but only genetically identified (Figure 1b). By quantifying the genomic divergence among these three groups, we can compare the relative effect of geographical isolation (allopatry) and natural selection (sympatry) on lineage divergence.
In this study, we aimed to (1) scrutinize the natural entity and divergence of TaiB, (2) clarify genomic and functional divergence among allopatric (Ryu-TaiB and Ryu-TaiA) and sympatric (TaiA-TaiB) groups, and (3) test the role of ghost introgression in sympatric divergence (TaiA-TaiB). Our framework provides a procedure for elucidating the evolutionary background of a sympatric cryptic lineage and highlights ghost introgression as a driver of speciation, particularly for species with low reproductive isolation and high extinction rate.

| MATERIAL S AND ME THODS
All analyses were performed using packages in R version 3.6.1 (R Core Team, 2019).

| Sampling and reference-based variant calling
The ddRAD raw data were retrieved from Chang et al. (2022).
Twenty-seven populations with 248 individuals, covering the whole distribution range, were examined. Raw reads were de-multiplexed following PCR duplicates (clone_filter) and phred score <30 reads removal by Stacks 2.53 (Catchen et al., 2013). The reference genome of C. panzhihuaensis (Liu, Wang, et al., 2022), a sister species with the same chromosome number (n = 11) as C. revoluta (Zonneveld, 2011), was used for reference-based mapping by Bowtie 2 (v2.3.4.1; Langmead & Salzberg, 2012). After coordinate sorting by picard (http://broad insti tute.github.io/picar d/, accessed on Feb 18, 2023), the bam files were processed in gstacks and populations programs in Stacks 2.53 for variant calling. Due to genetic admixture and lack of distinct groups, all individuals were treated as one population for the populations program (p = 1). Loci were retained if they were present in at least half of the individuals (r = .5) with minor allele frequency >5% (--min-maf = 0.05) and maximum observed heterozygosity <0.8 (--max-obs-het = 0.8). For further quality control by vcftools (Danecek et al., 2011), we set the total genotype depth and quality at <3 (--minDP = 3) and 10 (--minGQ = 10), respectively, as missing. Individuals having a missing rate >35% and loci with a missing rate >60% were removed. For more accurate island genetic diversity estimations, the island-specific missing SNPs were also removed.
Three datasets were created for further analyses: (1) containing all SNPs for detection of group-wise diverged F ST and ghost introgression genetic variations; (2) containing neutral SNPs for quantifying population genetic diversity, population structure, and divergence for the sympatric and allopatric groups; and (3) a dataset containing all SNPs of the Taiwan populations, TaiA and TaiB, for testing the evolution of the cryptic lineage after filtering missing SNP unique to TaiA and TaiB ( Table 1). Because of the admixed genetic structure, the genetic outliers were detected and removed by PCAdapt, which detects outliers without predefined grouping (Luu et al., 2017). The optimal K in the PCA scree plot was selected by Cattell's rule, and SNPs with −log p > 150 were defined as outliers.

| Population structure and divergence
Population structure was analyzed by sparse nonnegative matrix factorization (sNMF) and discriminant analysis of principal components (DAPC) in the LEA (Frichot et al., 2015) and adegenet packages (Jombart, 2008), respectively. In sNMF, the focus was on one global structure combining all individuals, and two separate structures, including Ryu + TaiA and TaiB. The cross-entropy of serial K from 1 to 15 was calculated for the best K selection with 1000 iterations. Q-matrices around the optimal K were plotted by pophelper (Francis, 2017). In DAPC, we only assessed the population structure of Ryu + TaiA because of the high divergence in TaiB. The best cluster was determined by the lowest BIC via the find.cluster function. The best number of PC for DA was retained by using optim.a.score function.
F I G U R E 1 (a) Phylogeny among Ryu + TaiA, and TaiB based on Chang et al. (2022). Colored circles match the classified groups based on geographic and genetic structure. (b) Sampling scheme. See Table S1 for detailed locality and population information. N and S groups are classified as two genetic groups in Chang et al. (2022). TaiA and TaiB refer to two sympatric lineages in CMR and Hongye populations in Taiwan. Ryu group includes Fujian, N, and S groups. CMR, coastal mountain range.
TaiB divergence was then estimated by an individual neighborjoining (NJ) tree and island pairwise F ST . The NJ tree was reconstructed by Nei's distance with 1000 bootstraps utilizing the aboot function in the poppr package (Kamvar et al., 2014). Pairwise F ST was calculated by arlsumstat v3.5.2 (Excoffier & Lischer, 2010).
We compared the divergence of TaiB with Ryu + TaiA populations.
Hierarchical clustering of F ST was also performed via the hclust function from stats after transforming by Euclidean distances.

| Population genetic diversity
Island-based and group-based (Ryu, TaiA, and TaiB) genetic diversity was calculated by a population program in Stacks 2.53. For estimating observed and expected heterozygosity (H obs and H exp ), Wright's inbreeding coefficient (F IS ), and private allele number (P N ) were considered. The private allele frequencies (P F ) of the cryptic, sympatric, and allopatric lineages were also compared.

| Group-wise divergence and functional annotation
In the dataset containing all SNPs (Dataset 1), three kinds of pairwise comparisons were calculated for detecting diverged SNPs between groups: (1) sympatric group with TaiA and TaiB, (2) allopatric group with TaiA and Ryu, and (3) allopatric group with TaiB and Ryu. Significantly diverged SNPs were defined by p < .001 in the Fisher's exact test, computed via the populations program (Catchen et al., 2013). A Venn diagram of the three groups was constructed to illustrate the group-wise divergence patterns.
The corrected AMOVA (analysis of molecular variance) F ST was used to compare sympatric and allopatric divergences: (1) rev versus tai: divergence between C. revoluta and C. taitungensis from previous taxonomy (TaiB was not included because its cryptic identity was not identified before); (2) Tai versus Ryu: allopatrically driven divergence between Taiwan and Ryu group; (3) TaiB versus Ryu: allopatrically driven divergence of TaiB from Ryu group; and (4) TaiB unique: SNPs that differentiate TaiB from allopatric Ryu and sympatric TaiA. The significantly diverged SNPs were annotated with the closest functional genes within ±1 kbp using C. panzhihuaensis as a reference through BEDOPS (Neph et al., 2012) and samtools (Li et al., 2009). Subsequently, the gene sequences were compared with the nonredundant (nr) database by the blastx-fast method in

| Evolutionary origin of cryptic lineages
Because of the apparent divergence between TaiB and Ryu + TaiA, we tested their evolutionary dependency by examining the current genetic exchange and the ghost introgression effect by hierarchical model testing using approximate Bayesian computation (ABC) with summary statistics (SumStats) and the random forest method (Pudlo et al., 2016;Raynal et al., 2019). With the ABC model settings, gene flow can be a good indicator for discerning allopatric and sympatric speciation in a species with poor dispersibility, like C. revoluta, from a population genetics perspective (Gavrilets, 2003;Richards et al., 2019). Gene-flow models between TaiA and TaiB were first tested for their evolutionary dependency. However, biased sample sizes between Ryu and TaiB may cause false-negative detection of gene flow. Thus, we only used the sympatric TaiA, which was genetically admixed with Ryu, to compare with TaiB in SumStats (Dataset 3).
If the two lineages are without current gene flow, complete isolation (CI) and primary contact (PC) would be supported; in contrast, secondary contact (SC) and continuous gene flow (CG) would be preferred. If the two lineages diverged allopatrically, gene flow would stop with the onset of isolation. Thus, CI or SC would be verified; alternatively, PC or CG is favored under sympatric divergence. Subsequently, we tested the ghost lineage models under the best of the four scenarios above. Ghost introgression may confer genetic variations in two ways: (1) gene flow between a hypothetical extinct ("ghost") lineage and TaiB or TaiA (scenarios GhGeB or GhGeA), or (2) hybridization between the ghost lineage and TaiB (or TaiA) to produce the hybrid lineage TaiA (or TaiB; scenarios GhHyA or GhHyB). The best gene-flow scenario without ghost gene flow (scenario GhPC) was compared with these four ghost models. The detected ghost introgression should arise from either the extinct species or unsampled extant species, which is now geographically distant from the current distribution of C. revoluta. Twenty-two summary statistics (SumStats) regarding population divergence, expansion, and size (Method S1) were applied in model selection and parameter estimation in the abcrf package (Raynal et al., 2019).

Dataset
Aims SNPs Individuals In model selection, the overlap of observed and simulated SumStats was first checked by PCA, then the sufficiency of the reference table size was evaluated by the stability of prior error rate, and posterior probability and models were validated for appropriate distinctions. Lastly, the best model was determined by the most positive votes, the posterior probability, and the goodness-of-fit (GOF) test against the observed SumStats.
In parameter estimation, the random forest hyperparameters were tuned first to find the proper node size and the number of variables. Second, the sufficiency of the reference table size  Tables S2 and S3; framework details in Figure S1 and Method S1).

| Ghost genomic introgression signal detection
Genomic signatures of introgressed regions from ghosts tended to show higher linkage disequilibrium (LD) and divergence than extant nonintrogressed regions (Li & Wu, 2021;Ottenburghs, 2020). Several haplotype-based methods have been developed to statistically search for these signals (e.g., sliding window, haplotype structure with LD pattern, etc.) (Li & Wu, 2021), but they are unsuitable for genomic markers like ddRAD that are based on reduced representation.
Consequently, we applied their rationale to identify the candidate introgressed functional genomic regions received from the ghost.
Based on Skov et al. (2018), the lineage receiving a ghost introgression will contain a high density of private SNPs. Therefore, we compared the genomic distribution of private SNPs, P F , F IS , and pairwise F ST between the ghost-introgressed and nonintrogressed lineages. As expected, the introgressed regions with the ghost showed clustered private SNPs, strong inbreeding, and genetic divergence with nonintrogressed lineage. If the ghost-introgressed genes underwent adaptive introgression, then the selective sweep should decrease the heterozygosity with high P F . Lastly, the candidate SNPs were functionally annotated using the method described above.

| Reference-based variant calling
After of 17% (Table 1). For Datasets 2 and 3, 860 and 801 SNPs were retained, respectively (Table 1). No individuals were removed from any of the three datasets. The absence of a "missingness" structure ( Figure 2) and the intermediate missing rate among cryptic lineage individuals (8%-20% and 5%-32% for all individuals; Figure 2) strongly supported the reliability of the datasets.
Moreover, all individuals of cryptic lineage (n = 10) came from four different ddRAD libraries, which ruled out the possibility of divergence caused by a batch effect.

| Distinct population structure and genetic divergence between Ryu + TaiA and TaiB
The number of clusters in sNMF was not determined in the global structure because K = 2-7 yielded similar cross-entropy states in Ryu + TaiA (mean = 0.036) signifies recent population expansion.
The small island populations were nearly panmictic, probably indicating they were sink populations that received genetic variations from multiple sources.

| Selective surpassed geographical process for cryptic lineage divergence
Comparisons between groups revealed an abundance of divergent SNPs between Ryu + TaiA and TaiB

| Ghost introgression and sympatric divergence
Reference  After hyperparametric tuning, a size of 50 nodes and seven variable splits generated the lowest prediction error. Table 3 and Figure S5 show the parameter estimation and posterior distribution for the

| Functional annotation implied stress resistance improvement and coastal adaptation of Ryu + TaiA by ghost introgression
Of 395 SNPs diverging in Ryu + TaiA and TaiB, 39 were unique to TaiB, and only four SNPs displayed ghost introgression signals (Table 4). These SNPs were densely localized on chromosome 2 F I G U R E 7 Hierarchical ABC and model selection results. Gray arrows indicate the gene-flow models as first tested, and the ghost lineage models were then tested based on the best gene-flow results. Red arrows highlight the differences between models. Parameters are described in Method S1 and results. Pos. prob, posterior probability of the best model. In gene-flow models: CI, complete isolation; CG, continuous gene flow; PC, primary contact; SC, secondary contact; in Ghost lineage models: GhGeA, ghost gene flow to A; GhGeB, ghost gene flow to B; GhHyA, ghost hybridization with A; GhHyB, ghost hybridization with B; GhPC, no ghost gene flow.

| DISCUSS ION
Previous research synonymized C. taitungensis as C. revoluta and reported the discovery of a cryptic lineage, TaiB, unique to Taiwan . In addition to the distinct genetic components

| Spatiotemporal background for ghost introgression
The best GhGeA model in ABC suggested ancient sympatric divergence of Ryu + TaiA and TaiB, occurring ca. 9 Mya in the late Miocene near the emergence of Taiwan at ca. 5-9 Mya (Sibuet & Hsu, 2004), assuming a generation time of 30 years . This result implied that the ancestral C. revoluta could have originated and diverged into Ryu + TaiA and TaiB in continental East Asia, where C. panzhihuaensis, the sister species of C. revoluta, was thriving Xiao & Moller, 2015).
Also, fossils resembling C. revoluta and C. panzhihuaensis have been discovered in this region (Su et al., 2014). During the late Miocene, there was an end to the global ocean circulation of warm tropical water, which resulted in a shift from a stable, warm climate to a cooler, more seasonal one (Steinthorsdottir et al., 2021). In East Asia, the changes were manifested latitudinally with cooler, drier weather in the northern regions and more seasonal, humid conditions in Indochina and southern East Asia (He et al., 2023), similar to current conditions where Cycas species have a high diversity (Lindstrom et al., 2009;Yessoufou et al., 2017). The crown age of the extant Cycas also dates back to the late Miocene, with a high diversification rate (Condamine et al., 2015;Liu, Wang, et al., 2022, but also see Liu et al., 2021 for Eocene diversification). These results indicated that the ancient Cycas diversity was high in the late

Miocene but may have retreated southward to southern East Asia
and Indochina with spatial contraction and demographic expansion or remained relatively stable in demography (Xiao & Moller, 2015).
The absence of fossil records after the late Miocene in northern East Asia also supports this perspective (Krassilov, 1978;Su et al., 2014;Yokoyama, 1911). Accordingly, as indicated by ABC, the ghost lineage may have been distributed around southern East Asia or Indochina when it went extinct or ceased the gene flow with C. revoluta during the middle Pleistocene.
Cycas species are known for their superior hybridization ability but are generally inferior in dispersal (Norstog & Nicholls, 1997;Whitelock, 2002). As a result, more condensed species richness and higher abundance of Cycas species diversity in southern East  following glaciation cycles (Kimura, 1996;Nakazawa & Bae, 2018).
Due to low dispersibility, ghost and ancestral C. revoluta would be contractions observed in many Cycas species (Feng et al., 2014;Gong et al., 2015;Liu et al., 2015;Zhan et al., 2011;Zheng et al., 2016). In summary, the GhGeA model in ABC suggests that the significant influence of ghost introgression on Ryu + TaiA and TaiB divergence since the late Miocene can be attributed to the paleoclimate-driven spatial dynamics of Cycas populations.

| Ghost introgression facilitates sympatric cryptic lineage divergence and ecological adaptation
Ghost introgression could facilitate adaptation and speciation by transferring advantageous genetic variations to its receivers.
The persistence of ghost introgression for ca. 8 million years has provided Ryu + TaiA with a rich source of adaptive genetic variations, potentially facilitating exaptation or local adaptation under novel environments. The colonization of C. revoluta in the Taiwan-Ryukyu Archipelago likely occurred through either a land bridge or when these islands were part of the Eurasian continental arc (Jiang et al., 2019;Li et al., 2017;Osozawa et al., 2012). Both eastward colonization routes exposed the populations to novel coastal environments from continent to island. Given the similarly large effective population size in Ryu + TaiA and TaiB, natural selection from the new coastal environment was likely to have been efficient for both lineages and accelerated their adaptive divergence.
Ghost introgression from the coastal-adapted lineage may have broadened the ecological resource spectrum of Ryu + TaiA and facilitated ecological adaptation from the ancestral inland to the coastal site (Germain et al., 2021). Spatial expansion may occur on these islands (Chiang et al., 2009). Although the nature of a reduced representation approach, such as ddRADseq, limits the detection of the ghost introgression signal, one functionally diverged candidate gene, HCHIB, related to ghost introgression was found.
Blister blight is a common worldwide problem for growers of C. revoluta (Faraz et al., 2020;Weber, 1944), with the symptom severity being positively associated with manganese (Mn) deficiency (Dehgan et al., 1994). The Mn deficiency is often linked to limestone areas (Kavvadias & Miller, 1999;Kowalenko & Ihnat, 2010), the primary habitats of the Ryukyu populations found in coastal regions.
Furthermore, coastal environments commonly experienced drought and heat stress, which were found to be associated with pleiotropic functions of HCHIB (Hayford et al., 2022), making populations with introgressed HCHIB better adapted to coastal habitats. Apart from HCHIB, other functionally diverged genes between Ryu + TaiA and TaiB have been identified across the genome, except for chromosome 6 (  Figure S4).
The absence of distinctive divergent SNPs between sympatric TaiA and TaiB (blue circles, Figure 6) also suggests that the divergence was not driven by the in situ emergence of newly derived SNPs. The significant drop in H obs (mean = 0.086) from H exp (mean = 0.14) among the three lineages raised concerns about severe inbreeding in TaiB (Figure 5a,b). Compared with the Ryukyu Archipelago and southern Kyushu, the habitat of C. revoluta in Taiwan was stable and more similar to its ancestral area in continental East Asia (Tsukada, 1966), serving as a refugium for preserving relic TaiB. However, for the Ryu + TaiA group, which would have an ecological shift from the continent ancestral populations, the Taiwan habitat could be a suboptimal niche for the marginal TaiA population, resulting in a similar inbreeding level but a lower private allele number with TaiB.

| Conservation of the invisibilities: cryptic lineage and ghost introgression
The recognition of cryptic lineages is important for accurate biodiversity estimation and effective conservation planning (Delic et al., 2017;Oury et al., 2021). However, caution must be exercised when dealing with recently derived cryptic lineages that lack reproductive isolation, as gene flow may hamper the detection of cryptic diversity and lead to taxonomic inflation (Chan et al., 2020(Chan et al., , 2021. Also, artificial genetic pollution from its congeners should be avoided. Contrarily, for the ancient divergent cryptic lineages with relic ancestral polymorphisms, such as TaiB of C. revoluta, protection of their refugial habitats is imperative. In Taiwan, the woody flora was comprised of paleo-and neo-endemics, forming a "mixed endemism" with high phylogenetic diversity and conservation value (López-Pujol et al., 2011;Wang et al., 2022). Correspondingly, the coexistence of relic TaiB and coastadapted TaiA in sympatry highlights the importance of conservation efforts in southern Taiwan (Cheng et al., 2005;Chiang et al., 2012;Huang et al., 2002;Kuo et al., 2009;Su et al., 2016;Wu et al., 2006).
Compared with TaiB, the flourishing Ryu + TaiA around the coast may have been introgressed by ghosts carrying the HCHIB resistance gene. Therefore, the study of ghost introgression can shed light on ancient ecological roles and inform conservation paleobiology (Dillon et al., 2022). Future research should integrate functional genomics, ancestral area reconstruction, and geohistorical and paleoecological data to enhance the beneficial effects of conservation efforts.

| CON CLUS IONS
The present study utilized reference-based ddRADseq to investigate the nature of sympatric cryptic lineages and evolution in C. revoluta. Our results revealed that ghost introgression was crucial in explaining sympatric divergence. Functional divergence of candidate genes from archaic gene flow implied adaptive introgression, facilitating ecological adaptation of Ryu + TaiA among islands from its continental ancestry. Moreover, the bottom-up functional annotations illuminated possible physiological and phenological differentiation, while top-down studies on the inconsistent morphological and ecological shift between Ryu + TaiA and TaiB can benefit taxonomic treatment and morphological evolution. Without careful morphological examination, early cryptic divergence in C. revoluta could have indicated morphological stasis following stabilizing selection, morphological convergence, or pseudo-cryptic divergence.
Overall, our findings shed light on the importance of ghost gene flow on extant species' adaptation and genomic divergence. The ghost would be the ancient cycad species or the unsampled species whose current distribution is distant from the focal species. While studies of genetic signatures and consequences of archaic introgression focus on humans due to their accessible archaic genome, we urge more attention to archaic introgression in nonmodel plants, particularly for those with less reproductive isolation than animals, to unravel the intricate web of life.

AUTH O R CO NTR I B UTI O N S
Jui-Tse Chang: Conceptualization (equal); formal analysis (equal); writing -original draft (equal); writing -review and editing (equal).

ACK N OWLED G M ENTS
We are grateful to Pei-Wei Sun for assistance with analyses. We